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O ! ABSTRACT 

u- 

We use the Galaxies-Intergalactic Medium Interaction Calculation (gimic) suite of 
i ^ i , cosmological hydrodynamical simulations to study the formation of stellar spheroids of 

Milky Way-mass disc galaxies. The simulations contain accurate treatments of metal- 
£NJ | dependent radiative cooling, star formation, supernova feedback, and chemodynamics, 

. and the large volumes that have been simulated yield an unprecedentedly large sam- 

' pie of pa 400 simulated ~ disc galaxies. The simulated galaxies are surrounded by 

(N ' low -mass, low-surface brightness stellar haloes that extend out to ~ 100 kpc and be- 

■ yond. The diffuse stellar distributions bear a remarkable resemblance to those observed 

around the Milky Way, M31 and other nearby galaxies, in terms of mass density, sur- 

C\J , face brightness, and metallicity profiles. We show that in situ star formation typically 

■ dominates the stellar spheroids by mass at radii of r < 30 kpc, whereas accretion of 

stars dominates at larger radii and this change in origin induces a change in slope of 
the surface brightness and metallicity profiles, which is also present in the observa- 
tional data. The system-to-system scatter in the in situ mass fractions of the spheroid, 
however, is large and spans over a factor of 4. Consequently, there is a large degree of 
scatter in the shape and normalisation of the spheroid density profile within r < 30 kpc 

^ ' (e.g., when fit by a spherical powerlaw profile the indices range from —2.6 to —3.4). 

We show that the in situ mass fraction of the spheroid is linked to the formation epoch 
of the system. Dynamically older systems have, on average, larger contributions from 
in situ star formation, although there is significant system-to-system scatter in this 
relationship. Thus, in situ star formation likely represents the solution to the long- 
standing failure of pure accretion-based models to reproduce the observed properties 
of the inner spheroid. 

Key words: Galaxy: evolution — Galaxy: formation — Galaxy: halo — galaxies: 
evolution — galaxies: formation — galaxies: haloes 



1 INTRODUCTION 

The stellar haloes of normal disc galaxies are expected to 
be excellent repositories of fossil evidence from the epoch 
of galaxy formation, even though these components con- 
tain only a small fraction of the total light. This is be- 
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cause the stellar haloes are widely believed to have been 
formed through the disruption of infalling satellite dwarf 
galaxies over the lifetime of the galaxy. Recent advances 
in observational techniques have allowed the detection of 
faint stellar haloes around external galaxies, such as M31 
jFerguson et al.ll2002l;lGuhathakurta et al-teOOotllrwin et ail 
20051: llbata et all 12007 ) and other nearby disc galaxies 
1 Mouhcine et all l2005al Tbl: Ide Jong et ail 120071 ; llbata et all 
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120091 ; iMouhcine et alll201Ch . Most stellar haloes have been 
shown to contain evidence of merger events in the form 
of tidal stellar streams. The number and physical prop- 
erties of these streams appear to be in gene ral agree- 
ment with the predictions o f th e ACDM model dBell et al.l 



20081: iGilbert et all l2009allbT: iMcConnachie et all 120091 ; 
Starkenburg et al.ll2009l ; iMartfnez-Delgado et al.ll2010h 



In spite of this, questions have been raised about the 
compatibility with the standard paradigm of the seemingly 
large differences in the properties of the stellar haloes and 
satellite systems of M31 and the Milky Way. These two sys- 
tems have similar masses but M31 appears to have expe- 
rienced a much more active merger history in the recent 
past, as suggested by its more disturbed stellar disc, larger 
bulge, younger and more metal-rich stellar populations, 
larger and mo re numerous tidal streams and surviving satel- 
lite galaxies jG uhathakurta et al. 20051 ; llrwin et al.l 2005; 
Brown et al 2006a h-jGilbert et alBood-lKalirai et al.lbood: 



McConnachie fc Irwi n 2006; Gilb ert et al. 120071 ; llbata et all 
20071 ; iKoch et al.ll2008l ; iKalirai et al.ll200g i. 

On the theoretical side, the detection of extended 
stellar components around external galaxies has reinvigo- 
rated the efforts of modelling stellar haloes in the con- 
text of ACDM cosmogony. Several theoretical studies sug- 
gest that the main properties of stellar haloes can be 



satellites dBullock & Johnston 20051; 


Font et all l2006allbllcl; 


De Lucia & Helmil 120081; iFont et all 


20081; IJohnston etail 


20081; Gilbert ct al. 2009bl; Cooper et al.l 20ld). These stud- 



ies have been carried out using a combination of high- 
resolution dark matter-only simulations and simple semi- 
analytic prescriptions for assigning light (i.e., stellar mass) 
to infalling sa tellites. Recent stu dies using such 'hybrid' 
methods (e.g., ICooper et ai1l201Gv i have taken advantage of 
the ultra-high resolution achieved by the lat est generation 
of da rk matte r only Milky Way simulations dSpringel et al.l 
120081 ; see also iDiemand. Kuhlen fc Madau Il2007h . which al- 
lowed them to explore the properties of even the faintest sub- 
structures (satellites and streams) surviving to the present 
day. 

However, even though a large fraction of the stel- 
lar hakQ of normal disc galaxies may be explained by 
the accretion scenario, there remain some stellar popula- 
tions in the Milky Way and other nearby galaxies that 
have not so far been accounted for in the context of 
accretion-only models. This has been known for a long time 
l|Eggen, Lvnden-Bell fc SandageHl962T ). but the evidence has 
become more convincing in the past few years. For ex- 
ample, the metal-rich halo stars in the Solar neighbour- 
hood have been shown t o have a bi-modal [a/Fe] pattern 
jNissen fc Schusterl [2010^ . In particular, the high [Fe/H] - 
high [a/Fe] population seems to be the product of very rapid 
chemical evolution that cannot be sustained in t he potential 
wells o f even the most massive accreted satellites |Font et al.l 
l2006ah . Instead, it has been suggested that these stars may 
have formed "in situ" in the halo (e.g.. IZolotov et al1l201Gi ) 
or they may have been stirred up from the disc by rep eated 
interactions with satellite galaxies l|Purcell et al.|[2010f ). Sim- 
ilarly, models that track only the accreted component of the 



1 By volume, but not by mass as we argue later in this paper. 



stellar halo (e.g.. IFont et al.ll200Sl ) have difficulties match- 
ing the intermediate age (~ 8 Gyr), metal-rich populations 
([Fe/ H] > -1) detected in the halo of M31 Csee lBrown et all 
l2007| y 

There may also be an indication of a spatial differ- 
entiation in the stellar halo of the Milky Way: analysis of 
main-sequence turnoff stars in the Sloan Digital Sky Survey 
(SDSS) sugge sts that the inner regions are more m etal-rich 
(e.g.. ICarollo et al.ll2007l , l20ld ; IdeJong et alj|20icl . but see 
ISesar et al.l 2011 ) and have a prograde motion with respect 
to the disc, while the outer halo is metal-poor and appears to 
have a small net retrograde rotation or no rot ation, depend 



ing on the adopted the local standard of r est jDeason et al 



2010; 



2011a) and/or the adop ted distance scale (|Schoenrich et al 



Beers et al.ll201lh . In addition, there is mounting evi- 



dence for a break in the slope of the density profile of the hal o 
beyond r ~ 25 kpc (Sesar ct al. 2 01 it iDeason et al.ll2011bT ). 
Furthermore, recent studies s uggest that the inner regions 
of th e halo are flattened ( e.g., lJurfc et "al1l2008l ; ISesar et"al] 
l201ll ; lDeason et alj[2011bf ). whic h is similar to what is seen 
for th e inner regions of M31 (e.g.. |Pritchet fc van den Berghl 
1 1994 ). These results suggest that the stellar haloes of galax- 
ies may have a dual nature, perh aps along the l ines o f 
the scenario proposed originally by ISearle fc Zinnl (|l97Sf ): 
while the outer halo formed primarily by accretion/mergers, 
some of the inner halo has formed through dissipational 
collapse. However, in the current paradigm, the dissipative 
comp onent would not have formed through monolithic col- 
lapse (|Eggen. Lvnden-Bell fc Sandagc 1962), but rather as 
a result of gas-rich mergers at early times. Gas-dynamical 
simulations of disc galaxies formed in a hierarchical sce- 
nario produce stellar systems with a significant in situ dis- 
sipative component (lAbadi et aL 20061; Brook et alj |2004| ; 
IZolotov et aLll2009l . l20ld : lOser et alj|20ich 7 It remains to be 
determined what the relative contributions of in situ and 
accreted stars are as a function of radius and system mass 
and, importantly, what the system-to-system scatter is in 
these contributions. The scatter in the in situ and accreted 
components may be at the heart of the reported large differ- 
ences in the properties of the extended stellar distributions 
of the Milky Way and M31. 

Another important aspect that we address here, which 
may be related to the in situ vs. accretion question, is the 
origin of large-scale metallicity gradients in stellar haloes. 
The emergence of metallicity gradients is a well-known con- 
sequence of dissipative collapse/mergers. The evidence in 
favour of metallicity gradients in local disc galaxies is mount- 
ing. For example, in the Milky Way, an underlying metal- 
licity gradient in the halo is suggested by the net differ- 
ence i n the mean metallicity between the inner and outer 
halo dEggen. Lvnden-Bell fc Sandagel 19621; Searle fc Zinnl 



19781; ICarollo et al.ll2007l . 120101 : Ide Jong et al.ll2010l . but see 



Sesar et al.ll201ll ). The halo of M31 also shows evidence for 



a gradient, where the inferred metallicity may drop by as 
much as 1 dex from the inner regions of the galaxy out to 
~ 100 kpc (see Fig. [5] below). NGC 891 represents another 
local Milky Way-analog for which there is evi dence of a stel- 
lar metallicity gradient (see llbata et al.l 2009; although note 
that metallicity measurements have only been made out to 
15 kpc for this system). 

Models that form haloes by non-dissipative accretion 
and mergers have so far not been able to produce significant 
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meta l licity gradients jFont et all l2006al ; |Pe Lucia fe Helmil 
120081 ; ICooper et all l2010l l. This suggests that either these 
models do not capture all the relevant physical processes 
(e.g., in situ star formation or ejection of disc stars by inter- 
action with satellites) or that they do not model the entire 
range of possible merger histories for Milky Way-mass galax- 
ies in the cosmological context. In this study we use cos- 
mological hydrodynamical simulations with a sophisticated 
treatment of chemodynamics to show that metallicity gra- 
dients are expected to be ubiquitous in the haloes of Milky 
Way-mass disc galaxies. The prominence of these gradients 
is shown to be tied to the fraction of the stellar mass that 
is formed in situ, which is itself determined by the merger 
histories of the galaxies. 

The paper is organised as follows. In Section [2] we de- 
scribe the simulations and the selection of the sample Milky 
Way-mass disc galaxies. In Section [3] we compute the av- 
erage properties of stellar haloes of these galaxies, such as 
their density profiles, metallicity gradients (both [Fe/H] and 
[a/Fe]) and metallicity distribution functions (MDFs). In 
Section|4]we quantify the contribution to the stellar spheroid 
from in situ star formation and accretion/disruption of satel- 
lites. Finally, we summarize and discuss our main findings 
in Section [5] 



2 SIMULATIONS 

We use the Galaxies-Intergalactic Medium Interaction Cal- 
culation (gimic) suite of cosmological hydrodynamical simu- 
lations, which were carrie d out by the Virgo C onsortium and 
are descr ibed in detail by Crain et al.l |2009l ) , hereafter C09 
(see also ICrain et al.l l 201ol . hereafter CIO). We will there- 
fore present only a brief summary of the simulations here, 
focusing on the details most relevant for the present study. 
This suite of simulations is ideal for the present study for 
several reasons: (i) it simulates large volumes of the universe 
that contain numerous Milky Way-mass disc galaxies, allow- 
ing us to robustly quantify both the mean trends and the 
scatter; (ii) it incorporates accurate physical prescriptions 
for metal-dependent radiative cooling, star formation, mass 
and energy feedback from Type la and Type II supernovae, 
as well as enrichment due to stellar evolution (AGB stars), 
which results in the formation of realistic disc galaxies; and 
(iii) the relatively high numerical resolution is adequate for 
resolving the stellar haloes of individual Milky Way-mass 
galaxies and their most massive satellite galaxies (i.e, the 
main contributors to the accreted component of the stellar 
halo) . 

The suite consists of a set of hydrodynamical re- 
simulations of five nearly spherical regions (~ 20h~ 1 
Mpc in radius) extrac ted from the Millennium Simulation 
jSpringel et al.l l2005aT ). The regions were selected to have 
overdensities at z = 1.5 that represent (+2, +1, 0, —1, — 2)cr, 
where a is the root-mean-square deviation from the mean 
on this spatial scale. The 5 spheres are therefore environ- 
mentally diverse in terms of the large-scale structure that is 
present within them. For the purposes of the present study, 
however, we will only select systems with total 'main halo' 
(i.e., the dominant subhalo in a friends-of-friends, hereafter 
FoF, group) masses similar to that of the Milky Way, ir- 
respective of which GIMIC volume it is in (i.e., irrespective 



of the large-scale environment). C09 found that the prop- 
erties of systems of fixed main halo mass do not depend 
significantly on the large-scale environment (see, e.g., Fig. 
8 of that study). We plan to study the properties of such 
galaxies as a function of their local environment in a future 
study. 

The cosmological parameters adopted for GIMIC are the 
same as those assumed in the Millennium Simulation and 
correspond to a ACDM model with Q. m — 0.25, Qa = 0.75, 
Qb = 0.045, as = 0.9 (where as is the rms amplitude of 
linearly evolved mass fluctuations on a scale of 8/i -1 Mpc 
at z = 0), H = 100/i km s~ x Mpc -1 , h = 0.73, n a = 1 
(where n s is the spectral index of the primordial power 
spectrum) . The value of as is approximately 2-sigma higher 
than has been inferr ed from the most recent CMB data 
l|Komatsu et al.ll2009l ). which will affect the abundance of 
Milky Way-mass systems somewhat, but should not signifi- 
cantly affect their individual properties. 

The simulations were evolved to z = using the 
TreeP M-SPH code gadget-3 (last described in ISpringel I 
2005). The gadget-3 code has been substantially modi- 
fied to incorporate new baryonic physics. Radiative cooling 
rates for the gas are computed on an element-by-element ba- 
sis by interpolat ing within pre-comp uted tables generated 
with CLOUDY (|Ferland et alj 1 1998 ) that contain cooling 
rates as a function of density, temperature, and redshift and 
that account for the presence of th e cosmic microwave back - 
ground and photoionisation from a lHaardt fc Madaul (l200lh 
ionising UV/X-Ray background fsee lWiersma et alj|2009af ). 
This background is switched on at z = 9 in the simulation, 
where it 'reionises' the entire simulation volume. Star for- 
mation is_tnickedjri_j4i^ the prescrip- 
tion of ISchave fc Dalla Vecchial ([2008) . Gas with densities 
exceeding the critical density for the onset of the thermo- 
gravitational i nstability is e xpected to be multiphase and 
to form stars |Schave 1120041 ). Because the simulations lack 
both the physics and the resolution to model the cold inter- 
stellar gas phase, an effective equation of state (EOS) is im- 
posed with pressure P cx /9 4//3 for d ensitiefl uh > n* where 
0.1 cm -3 . As described in ISchave fc Dalla Vecchial 



(2008), gas on the effective EOS is allowed to form stars 
at a pressure-dependent rate that reprod uces the observed 
Kennicutt-Schmidt law ( Ke nnicuttl 1 19981 ) by construction. 
The timed release of individual elements by both massive 
(Type II SNe and stellar winds) and intermediate mass stars 
(Type la SNe and asymptotic giant branch stars) is i ncluded 
following the prescription of IWiersma et ah! l|2009bl ). A set 
of 11 individual elements are followed in these simulations 
(H, He, C, Ca, N, O, Ne, Mg, S, Si, Fe), which represent all 
the important species for computing radiative cooling rates. 

Feedback from supernovae is incorporated usi ng the ki- 
netic wind model of lDalla Vecchia fc Schavel (|2008l ) with the 
initial wind velocity, v w , set to 600 km/s and the mass- 
loading parameter (i.e., the ratio of the mass of gas given 
a velocity kick to that turned into newly formed star parti- 



2 Gas particles are only placed on the EOS if their temperature 
is below 10 5 K when they cross the density threshold and if their 
density exceeds 57.7 times the cosmic mean. These criteria pre- 
vent star formation in intraclust er gas and in intergalactic gas at 
very high redshift, respectively (Sc have fc Dalla Vec chia 2008). 
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cles), 77, set to 4. This corresponds to using approximately 
80% of the tot al energy available from supernovae for a 
IChabrierl (2003) IMF, which is assumed in the simulation. 
This choice of parameters results in a good match to the 
peak of the star formation rate history of the universe (C09; 
see also ISchave et afll20ld ) and reproduces a number of X- 
ray/optical scaling relations for normal disc galaxies (CIO). 

The GIMIC simulations have been run at three levels of 
numerical resolution, 'low', 'intermediate', and 'high', to test 
for numerical convergence. The low-resolution simulations 
have the same mass resolution as the Millennium Simula- 
tion (which, however, contains only dark matter) while the 
intermediate- and high-resolution simulations have 8 and 64 
times better mass resolution, respectively. In the case of the 
high-resolution simulations, only the —2a volume was run to 
z = 0, owing to the computational expense of running such 
large volumes at this resolution. For the purposes of the 
present study, we therefore follow the approach of CIO and 
adopt the intermediate-resolution simulations in the main 
analysis and reserve the high-resolution — 2a simulation to 
test the numerical convergence of our results. 

The intermediate-resolution runs therefore have a dark 
matter particle mass moM — 5.30 x lO 7 /^ 1 M and an 
initial gas particle mass of m gas ~ 1.16 x 10 7 /i _1 Mq, im- 
plying that it is possible to resolve satellites with total stel- 
lar masses similar to those of the classical dwarf galaxies 
around the Milky Way (with Af» ~ 10 9 " 10 M Q ) with several 
hundred up to ~ 1000 particles. This is important for the 
present study, as it is currently believed that the disruption 
of such massive satellites is the primary contributor by mass 
to th e accreted stellar halo (|Font et al.ll2006al ; ICooper et al.l 
2010). Without a sufficiently large number of particles, the 
internal structure of the satellites will not be adequately re- 
solved which can lead to overefficient tidal stripping of the 
satellites. In addition, relatively high resolution is required 
to model star formation in th e satellites correctly . For e xam- 
ple, iBrooks et all |2007l ) and IChristensen et all l|2010t ) find 
they require at least several thousand to ~ 10 4 particles 
before convergent star formation histories are obtained for 
isolated dwarf galaxies simulated with the SPH code GASO- 
LINE. This is similar to the number of particles with which 
we resolve the most massive satellites galaxies in the high- 
resolution GIMIC run. In the Appendix, we present a con- 
vergence study for our simulations, comparing the results of 
the intermediate- and high-resolution GIMIC simulations. We 
show that the present-day luminosity function of satellites, 
as well as the mass and metal distribution of the stellar halo, 
in our default intermediate-resolution runs are robust to an 
increase in the mass resolution of a factor of 8. In Table [1] 
we present a summary of the mass and force resolution of 
the simulations used in this study. 

Note that for all these simulations the remainder of the 
(500/1" 1 Mpc) 3 Millennium volume is also simulated, but 
with dark matter only and at much lower resolution. This en- 
sures that the influence of the surrounding large-scale struc- 
ture is accurately accounted for. 

2.1 Selection of Milky Way- mass systems 

Below we describe the selection criteria for our sample of 
Milky Way-mass systems. 

First we define what we mean by a 'Milky Way- 



Table 1. Resolution of the GIMIC runs. JVdm and jV* are the 
median number of dark matter and star particles within r2oo °f 
the simulated galaxies, ttidm and m gaa are the dark matter and 
(initial) gas particle masses. e so f t is the Plummer equivalent force 
resolution in physical space at z ^ 3. 



Run 


JVdm 


iV* 


m DM 


"^gas 


e soft 




(10 4 ) 


(10 4 ) 


(M Q /h) 


(M e /h) 


(kpc/ft) 


Int-res 


1.6 


0.59 


5.30 x 10 7 


1.17 x 10 7 


1.0 


Hi-res 


9.1 


3.7 


6.63 x 10 6 


1.46 x 10 6 


0.5 



mass' system. This is a system with a present-day total 
(gas+stars+dark matter) mass within V2oo (i-e., the radius 
which encloses a mean density of 200 times the critical den- 
sity of the universe at z = 0) of 7 x 10 U M Q < Af 20 o < 
3 x 10 12 Mq. For reference, this range roughly spans the mass 
estimates in the l iterature for both the Milky Way and M31 
[|Kochaneklll996l: lEvans fc Wilkinsonl l200d: iBattaglia et all 
120051 ; iKaranchentsev fc Kashibadzell2006l ; iGuo et al.ll2QloK 
although we are not attempting to reproduce the properties 
of either of these systems in detail. Systems and their sub- 
structures are id entified in the simu lations using the Sub- 
find algorithm o f Dolag et al. I(l2009h. that extends the stan- 
dard implementation of ISpringel et all l)200ll ) by including 
baryonic particles in the identification of self-bound sub- 
structures. For each FoF system that is identified, a spher- 
ical overdensity (SO) mass with A = 200 is computed 
(i.e., M2oo)- The gas+stars+dark matter associated with 
the most massive subhalo of a FoF system are considered 
to belong to the 'galaxy', while the gas+stars+dark matter 
associated with the other substructures are classified as be- 
longing to satellite galaxies. As the current study is focused 
on examining the extended stellar distributions of normal 
disc galaxies, we exclude all bound substructures (i.e., satel- 
lites) from the analyses presented in this paper, with the 
exception of Fig. [2] (described below) . 

Galaxies identified in this way are then classified mor- 
phologically as disc- or spheroid-dominated, based on their 
dynamics. In particular, we use the ratio of disc stellar mass 
to total stellar mass (D/T) within 20 kpc computed by C10. 
We select disc galaxies with D/T > 0.3, although very simi- 
lar results are obtained for other cuts in D/T (e.g., we have 
also tried 0.2 and 0.4 and none of the main results of our 
study are changed). C10 decomposed the spheroid compo- 
nent from the disc by first computing the angular momen- 
tum vector of all the stars within 20 kpc, calculating the 
mass of stars that are counter-rotating, and doubling it. This 
procedure therefore assumes the spheroid has no net angu- 
lar momentum. The remaining mass of stars is assigned to 
the disc and D/T is computed as the disc mass divided by 
the sum of the disc and spheroid masses. C10 plotted the 
D/T distribution in their Fig. 1 for a sample of galaxies in 
the intermediate resolution GIMIC runs with a very similar 
mass range to the one adopted in the present study. They 
found that approximately half of all the simulated Milky 
Way- mass galaxies are disc-dominated (D/T > 0.5), while 
approxiately 20% lack a significant spheroidal component 
D/T > 0.75 (i.e., are effectively "bulgeless", Sales et al. in 
prep). 

The method of C10 for decomposing the disc and 
spheroid does not assign individual star particles to either 
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Figure 1. Surface mass density maps of all the stars (left) and gas (right) in one of the sample Milky Way-mass disc galaxies at z = 
in the high-resolution simulation. The panel is 100 kpc on a side. 



component, it simply computes the total mass in two com- 
ponen ts. We therefore follow the approach of lAbadi et al.l 
(|2003h in assigning particles to the disc and spheroid. In 
particular, we align the total angular momentum of the stel- 
lar disc with the z-axis, calculate the angular momentum of 
each star particle in the x-y plane, J z , and compare it with 
the angular momentum, J c irc, of a star particle on a co- 
rotating circular orbit with the same energy. Disc particles 
are se lected using a cut of J z /J c irc > 0.8 as in lZolotov et al.l 

(2009) . Visually, this cut is quite successful in isolating the 
disc, but in a small number of cases we found that particles 
well above/below the disc mid-plane were also assigned to 
the disc. For this reason, we also apply a spatial cut such that 
disc particles cannot exist more than two softening lengths 
above or below the disc mid-plane. Note that we have also 
tried assigning disc particles based on the fracti on of stellar 
kineti c energy in ordered rotation, as defined in lSales et al.l 

(2010) , but found no significant differences from the method 
outlined above. 

One can al s o com pute the D/T ratio using the method 
of lAbadi et al.l (I2003T) with the angular momentum cut of 
IZolotov et al. l|2009l ). Comparison of the D/T values com- 
puted this way to that of the method of C10 shows a strong 
correlation, in the sense that both methods agree on which 
galaxies are most 'disky'. In general, though, the D/T ratio 
computed using the method of C10 is m 0.2 larger on average 
than that computed using the method of lAbadi et al.l (120031 ') 
with the angular momentum cut of IZolotov et al.l ([2009) . 
The implication of this is that the spheroid does have net 
angular momentum in reality (or else that the stellar disc 
has a much broader range of an gular momenta than that 
advocated bv lZolotov et al.ll2009l ). As we have noted above, 
however, our main conclusions are insensitive to the exact 
cut in D/T adopted for our galaxy sample selection. 



Unlike lZolotov et all |2009l ). we have elected not to fur- 
ther decompose the spheroid into bulge and halo compo- 
nents. While separating the disc from the spheroid is rela- 
tively straightforward to do for both the observations and 
the simulations, separating the bulge from the stellar halo is 
more challenging, as there is no sharp distinction in the dy- 
namics or the spatial distribution of the two components. A 
variety of techniques have been used in theoretical and ob- 
servational studies to distinguish these components, which 
can complicate the comparison of any single component be- 
tween different studies (at least in the inner regions where 
the bulge and halo overlap spatially). Our approach is to 
avoid this complication by simply reporting on the total 
bulge+halo component, which can readily be deduced from 
observations. 

Lastly, we point out that we have not imposed any ex- 
plicit constraints on the merger histories of the simulated 
galaxies. This is in contrast with many of the previous 
cosmological modelling studies of Milky Way-type galax- 
ies, which usually select their simulation candidates from 
systems that have n ot undergone any major mergers in 
the recent past (e.g ., Diemand et al.l 120051 ; lOkamoto et al.l 
l2005l;ICooper et alJ|201Ch . This condition is meant to ensure 
that these systems develop stellar discs that survive until 
the present day and t hat resemble the Milky Way's disc 
l|T6th fc Ostrikerlll992l ). Note that these previous studies 
were based on hybrid methods applied to dark matter-only 
simulations and therefore did not have the benefit of be- 
ing able to check which histories actually give rise to discs 
at the present day. Importantly, more than 70% of Milky 
Way-mass haloes have experienced major mergers (defined 
as Msat/Mdisc ^ 1/3, where M sa t is the total mass of the 
satellite and Mdi 3C is the stell a r mass of the disc) in the las t 
10 Gyr l|Stewart et al.l |2008| ; iBovlan-Kolchin et all I2010T I, 
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Table 2. Global properties of the simulated disc galaxies. Presented is the median of the property in 3 different total mass (M200) bins, 
with the error bars representing the 5 th and 95 th percentiles. My is total V-band magnitude, M„(< r20o) is the total stellar mass within 
T200, UrotCR©) is the rotation velocity of the stellar disc at the Solar radius, [Fe/H] r< 30k pc is the mean stellar metallicity within 30 kpc, 
[Fe/H] 7 . > 3ok pc is the mean stellar metallicity beyond 30 kpc, and n g!i i,bin is the number of simulated galaxies in the total mass bin. 
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implying that the majority of merger histories may be ex- 
cluded from studies that select only quiescent histories. The 
constraints applied to previous hybrid studies appear to be 
overly conservative, given that the majority of the Milky 
Way-mass systems in GIMIC are disc-dominated systems (see 
C10). The resiliency of galactic discs to mergers likely stems 
from the fact that the galaxies contain a significant dissipa- 
tional gaseous component (th at was not taken into account 
in, e.g., Toth &: Ostr ikcr 1992 ), which is able t o maintain (or 
re-establish) the disc (see lHopkins et al"1l2009l ). In any case, 
for the present study, which is based on hydrodynamic sim- 
ulations of large volumes, there is no need to place explicit 
constraints on the merger histories since we know which 
galaxies have prominent discs at z = 0. 

Applying the system mass and morphology criteria de- 
scribed above yields a total sample of 412 systems in the 5 
intermediate-resolution GIMIC volumes. This greatly exceeds 
the sample size of all previous studies of this type (i.e., those 
based on hydrodynamic re-simulations or hybrid methods), 
al beit at lower resolut ion than some previous studies (such 
as IZolotov et al]|2009u . The large number of systems in our 
sample allows us to robustly quantify the system-to-system 
scatter present in the properties of stellar halo es. Our work 
is the refore complementa r y to s tudies such as lAbadi et al.l 
|2006t ) and lZolotov et all (2009), which are based on high- 
resolution simulations of a small number of systems. (But 
note we demonstrate in the Appendix that the properties of 
our simulated stellar haloes are robust to increased numeri- 
cal resolution.) 

As an example, Fig. [T] shows the stellar and cold gas 
mass distributions of a typical simulated disc galaxy. An 
extended, flattened stellar distribution is clearly visible out 
to large radii. 



3 PROPERTIES OF STELLAR HALOES 

Before analysing the detailed properties of the stellar haloes 
of the simulated disc galaxies, we briefly discuss the overall 
global properties of the sample galaxies. The motivation for 
this is simple: if the simulated galaxies were to look wholly 
unlike observed normal disc galaxies, then that would limit 
what can be concluded about the formation of their stel- 
lar haloes. As we show in Section 3.1, the simulated galax- 
ies have global properties (i.e., total luminosities and stel- 
lar masses, rotation velocities, and satellite luminosity func- 
tions) which are comparable to those estimated for the Milky 



Way and M31. The reader who is interested only in the prop- 
erties of the stellar haloes of the simulated galaxies may wish 
to skip ahead to Section 3.2. 

3.1 Global properties 

In Table 1 we present the typical properties of the simulated 
disc galaxies, split into several halo mass (M200) bins. In par- 
ticular, we list total V-band magnitude (My, in AB mags.), 
total stellar mas^f] within r2oo [M*(< raoo)], disc rotation 
speed at the solar radius [v ro t(RQ), where Rq is assumed to 
be 8.5 kpc] and mean stellar metallicity within and beyond 
a galacto-centric distance of 30 kpc ([Fe/H] r< 30k pc )- (The 
choice of 30 kpc is motivated below, where we show how the 
physical nature of the bulge+halo component changes with 
radius.) The values listed in Table 1 are medians for galaxies 
in the mass bin, while the error bars represent the 5 th and 
95* percentiles. We also list the number of galaxies in each 
mass bin (% a i,bln)- 

The V-band luminosity of star particles in the simula- 
tions is calculated as follows. Each star particle is treated as 
a simple stellar population (SSP). The ages and metallici- 
ties of star particles are used to compute the spectral energy 
distri bution (SEP) by interpolat ing over the Galaxev mod- 
els of iBruzual fc Chariot! (|2003h . The V-band luminosity is 
obtained by integrating the product of the SED with the 
V-band filter transmission function. We ignore the effects of 
dust attenuation in this calculation. 

As already noted, the GIMIC simulations successfully re- 
produce the peak of the observed star formation rate density 
history of the Universe, as well as the X-ray-optical scaling 
properties of normal disc galaxies (C10). Another notable 
success of the simulations is that they produce large num- 
bers of disc galaxies (see Fig. 1 of C10). Realistic spiral discs 
have been notoriously difficult to achieve in the past with 
cosmological hydrodynamic simulations. We will present a 
detailed analysis of the galaxy discs in a future study (Mc- 
Carthy et al., in prep.). 

The GIMIC simulations, however, do not reproduce 
the bright end of the global galaxy luminosity function, 

3 When referring to stellar masses throughout the paper we are 
always referring to the current stellar mass, which is less than 
the initial stellar mass due to stellar mass loss. Typically, about 
45% of the initial mass is lost due to stellar evolution, which is 
dictated by the Chabrier IMF adopted in the simulations and the 
generally old stellar populations we are studying here. 
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Figure 2. The mean satellite luminosity function of simulated 
galaxies in the intermediate-resolution GIMIC simulation. The 
black dashed curve represents the mean luminosity function for 
all the simulated galaxies. The thick black solid curve represents 
the mean luminosity function for simulated galaxies in the mass 
range 1.8 X 10 12 < M 200 /Mq < 3.2 X 10 12 (128 disc galaxies meet 
this criteria). The hatched region represents the Poisson error 
distribution for this luminosity function. The solid, red and blue 
histograms show the satellite luminosity functions for the Milky 
Way and M3 1, respectively. The Milky Way data is taken from 
iMated l|l998l l and the M31 data is taken fromlMcConnachi e et all 
(2009). Overall, massive galaxies in our sample produce luminos- 
ity functions that are quite similar to those observed for the Milky 
Way and M31, suggesting that the accreted component of the halo 
in the simulations should have a realistic mass fraction. 



Mv < —22 (see C09). This is perhaps not surpris- 
ing, as the simulations lack feedback from supermassive 
black holes, which is thought to be crucial for regulat- 
ing star formation in the most massive systems (e.g. , 
Springel et al.ll2005bl; [Bower et al.|[200r3 ; ICroton et al]|2006l: 



Booth fc Schavell2009r r The match to the observed abun- 
dance of ~ L, galaxies is, however, reasonable. One should 
also recall that the luminosity function depends not only 
on the luminosity as a function of halo mass (which we 
quote below), but also on the number density of dark matter 
haloes, which depends on the adopted cosmological param- 
eters. This implies that one can obtain correct properties of 
individual galaxies but still not reproduce the observed lu- 
minosity function, if the wrong cosmological parameters are 
adopted. Therefore, a better approach to assessing whether 
the simulated galaxies are realistic, is to compare the mass- 
to- light ratio (or total luminosity) as a function of halo mass 
with the observations, which we do immediately below. 

For comparison, M31 has an absolute visual magnitude 
of M v w -21.0, M* ~ 5.6 x 10 10 A/ o , and v Iot « 250 km/s. 
The absolute magnitude was calculated using the GALEX 
apparent V magnitude of iGil de Paz et ail |2007r ) and as- 
suming a distance to M31 of 785 kp c. The total stellar mass 
was calculated using B-V and V of iGil de Paz et all (|2007i ) 



to infe r the stellar M/L using the relation of lBell fc de Jong] 
|200ll ) and scaling the resul t by 0.7 to c o nvert from 
the Salpeter IMF assumed by iBell fc de Jond ([200 If ) to a 
Chabrier IMF, which is adopted in the GIMIC simulations. 
The rotation velocity, which has been corrected for incli- 
nation, was obtained from the HyperLeda data basfl The 
Milky Way h as a total stellar mass of M , « 5.5 x 10 10 M (7) 
l|Flvnn et al.ll200 6fl and v rot = 250 km/s (|Reid et aUl2009h . 
These properties of M31 and the Milky Way are quite typ- 
ical of the simulated galaxies, particularly those with halo 
masses of A'hoo < 2 x 10 12 Mq. We will comment in detail 
on the metallicity of the stellar halo of both observed and 
simulated systems below. 

How do the D/T ratios for M31 and the Milky Way 
compare with that of our simulated galaxies? The D/T ra- 
tio for M31 and th e Mi lky Way have been estim ated by 
iGeehan et ail l|2006h and lDehnen fc Binnevl (|l998l ), respec- 
tively, by modelling the surface brighrtness distribution of 
the g alaxies and ass u ming a particular stellar mass-to-light 



ratio. 



Geehan et al 



iDehnen fc Binnevl (fl998i rfind D/T 



(2006) find D/T 



=a 0.69 for M31 and 
0.73 - 0.92 (the ex- 
act value depending on the surface brightness modelling as- 
sumptions). These values fall within the simulated D/T dis- 
tribution when the dynamical method of C10 is used but 
they are la rger than any of val ues of D/T estimated using the 
method of lAbadi et al 1 (I2003T I with the angular momentum 
cut of IZolotov et alj (120091 ) . This begs the question, which 
is the appropriate comparison to make? The answer is likely 
to be neither. As pointed out recently bv lScannapieco et al.l 
(2010), there can be large systematic differences in the D/T 
ratios inferred through kinematic and 'morphological' (i.e., 
based on surface brightness modelling) methods. In partic- 
ular, those authors found by analysing simulated galaxies in 
the same way as observers do (i.e., morphologically, includ- 
ing dust attenuation) they typically recovered much higher 
D/T ratios than that recovered with the standard dynam- 
ical decomposition applied by simulators. Thus, to make a 
more meaningful comparison between our simulated galax- 
ies and M31 and the Milky Way would require analysing our 
simulations in the same way as has been done for these two 
galaxies. This is clearly a worthwhile exercise but is beyond 
the scope of the present study. 

Since the present work concerns itself with the prop- 
erties of the extended stellar distributions of disc galaxies 
(i.e., the spheroids), it is also important to check whether 
the simulations yield reasonable satellite luminosity func- 
tions and radial distributions, since a significant fraction of 
the spheroid (by volume at least) is thought to have been 
assembled through the tidal disruption of massive satellites. 
In Fig. [5]we compare the mean satellite luminosity functions 
(LFs) of simulated galaxies at z = (both for all galaxies in 
our sample and for a subset of high-mass simulated galax- 
ies) with those of the Milky Way and M31. For the simu- 
lated galaxies the LF is calculated for all satellites within 300 
kpc, similar to the outermost radius adopted for the obser- 
vational data. The LFs of the satellite systems of the Milky 



4 http: / /leda. u niv-ly onl.fr. 

5 iFlvnn et all d2006h directly measure the stellar M/L locally 
and apply it to the Galaxy to estimate the total stellar mass. The 
measured value of M/L is consistent with the predictions of a 
Chabrier IMF. 
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Way and M31 have a slightly higher normalisation than the 
mean LF derived from all the simulated galaxies. However, 
encouragingly, they are quite compatible with a large sub- 
set of massive galaxies in our sample, with masses in the 
range 1.8 x 10 12 < M 200 /Mq < 3.2 x 10 12 (128 disc galax- 
ies meet this criteria). The observationally inferred masses 
for t he MW and M31 are fully consistent with this range 
(e.g., iGuo et al.ll20ld ). Note that the flattening in the LF 
of the simulated satellites at My > —12.5 is due to the 
finite resolution of our simulations (see the Appendix for 
a convergence study ). In terms of the spatial distribution, 
iDeason et al.l ((201 ld ^) examined the cumulative radial distri- 
bution of the 10 brightest satellites around Milky Way-mass 
galaxies in the intermediate-resolution GIMIC simulations. 
They find that the Milky Way and M31 distributions fall 
within the system-to-system scatter of the simulated galax- 
ies. Taken together, these results imply that the simulated 
spheroids should have realistic contributions from the accre- 
tion and disruption of satellites. 

We now proceed to our analysis of the stellar spheroids 
of simulated disc galaxies. 



3.2 Structural properties 

Below we present the radial distributions of stellar mass den- 
sity, p*, and surface brightness, Ey, as well as the radial 
profiles of the metallicity [Fe/H] and [a/Fe] of our sample 
of simulated galaxies. 

3.2.1 Stellar mass density and surface brightness profiles 

Fig. [3] shows the median spherically-averaged stellar mass 
density (left panel) and azimuthally-averaged V-band sur- 
face brightness (right panel) radial profiles of the bulge+halo 
component (solid black curves), as well as the 1- and 2- 
sigma scatter about the median (shaded regions). The solid 
red curves show the median profiles for the total stellar 
(disc+bulge+halo) components. The dashed green curve 
represents a power-law in p, with index of —3.5 in the left 
panel, which corresponds to a power-law in surface bright- 
ness (in linear units) with index of —2.5, assuming the stel- 
lar mass-to-light ratio is independent of radius. Converting 
to mags, per arcsec 2 , implies a power-law with index 5.25, 
which is represented by the dashed green curve in the right 
panel. 

First, extended (out to r ~ 100 kpc and beyond) stel- 
lar distributions are a ubiquitous feature of the simulated 
galaxy populations. The radial distributions of mass den- 
sity and surface brightness have relatively simple forms, but 
as a comparison of the solid black and dashed green curves 
demonstrates, no single power-law distribution can match 
the profiles over all radii. A power-law with index of —3.5 
(in p„) provides an excellent fit to the outer regions, r > 30 
kpc, of the median bulge+halo profile, but the same power- 
law under-predicts the mass density/surface brightness at 
smaller radii. As we will discuss in Section 0] the change 
in the nature of the radial distributions at radii of r « 30 
kpc is driven by a change in the relative contributions of 
stars deposited by tidally disrupted satellites and stars that 
formed in situ. 

We can quantify the system-to-system scatter in the 



slope of the profile at large radii by fitting power-laws to the 
mass density distributions of all « 400 disc galaxies in our 
sample. We find a median slope of —3.51 with —3.05, —3.20, 
—3.89, and —4.16 as the 5th, 14th, 86th, and 95th percentiles 
of the distribution of slopes, respectively. Although it usually 
does not fit the profiles as well as a power-law, we have also 
fit a de Vaucouleurs profi le to the outer regions (as done, 
e.g., bv llbata et al.1 .2007 for the outer regions of M31). In 
this case we find a median effective radius, r e , of 17.80 kpc 
with 6.71, 10.05, 34.74, and 54.98 kpc as the 5th, 14th, 86th, 
and 95th percentiles, respectively. 

In the inner regions, r < 30 kpc, we find that a de 
Vaucouleurs profile usually fits the mass density distribution 
very well. Here we find a median r e of 1.83 kpc with 0.73, 
0.95, 4.90, and 7.54 kpc as the 5th, 14th, 86th, and 95th 
percentiles, respectively. A power-law fit to r < 30 kpc yields 
a median slope of —2.98 and —2.47, —2.58, —3.30, and —3.36 
as the 5th, 14th, 86th, and 95th percentiles, respectively. 
This is shallower than in the outer regions, but note that 
the mass density profile has a higher normalisation at small 
radii, so that a power-law with slope —3.5 extrapolated to 
the inner regions underpredicts the mass density there. 

In general, our results confirm those of several previous 
theoretical studies on the sha pe of the mass density pro - 
file at larger radii. For example, iBullock fc Johnston! l(2005h . 
who used a hybrid semi-analytic plus idealised N-body ap- 
proach to study the formation of the stellar halo, find that 
a power-law with index of —3.5 describes the shape of their 
haloes very well at radii of ~ 50 — 100 kpc. More recently, 
ICooper et all (|201fjh applied a similar approach to the very 
high-resolution Aquarius cosmological N-body simulations. 
These authors did not report on a power-law fit to their 
simulated profiles at large radii, but an examination of their 
Fig. 4 suggests a slightly steeper logarithmic slope of ~ —4. 
The shape of the bulge+ halo component at large radii is also 
similar to that found bv lAbadi et al.l l|2006l ). who performed 
cosmological SPH re-simulations of a small number of Milky 
Way-mass systems. In particular, these authors found log- 
arithmic slopes of —2.9 at 100 kpc and —3.5 at the virial 
radius of their simulated galaxies. 

While there is generally good agreement at large radii 
between our results and those of previous theoretical stud- 
ies, som e differences become appare nt at small radii. In par- 
ticular, IBullock fc Johnston! (|2005l l find that their profiles 
become increasingly flat at small radii, so that the loga- 
rithmic slope approaches —1 for r ~ 10 kpc, whereas our 
spheroid component shows no evidence of flattening at small 
radii - on the contrary, the profile actually steepens some- 
what at radii of r m 30 kpc before leveling off to a slope 
of sa —3 at r < 20 kpc. This dramatic flattening in the 
IBullock fc Johnston! (120051 ) model is almost certainly in part 
due to the neglect of in situ star formation in their hybrid 
N-body model, but another potentially important factor is 
the treatment of the galaxy potential well, which does not 
feel the back-reaction due to the orbiting sat ellites in their 
idealis ed N-bod y simulations. For e xample, ICooper et al.l 
(2010) (see also lDiemand et alj|2005h . who applied a simi- 
lar hybrid method to cosmological N-body simulations with 
'live' galaxy potential wells, find a relatively steep distribu- 
tion at small radii, with a logarithmic slope of ~ —3. This 
suggests that the treatment of the potential well is indeed 
relevant. Consistent with this hypothesis is that when we fit 



Stellar haloes of disc galaxies 9 




Figure 3. Left: Median spherically-averaged mass density profiles. The solid red curve is the median stellar mass density profile for 
all stars, whereas the solid black curve is the median stellar mass density profile for just the spheroidal (bulgc+halo) component. The 
finely hatched region encloses the 14 th and 86 th percentiles of the stellar bulge+halo mass density profile, while the sparsely hatched 
region encloses the 5 th and 95 th percentiles. The dashed green curve represents a power-law profile with p„ <x r -3 5 with a normalisation 
chosen to match the median stellar halo mass density profile at large rad ii. The dot-dashed blue curve r epresents the median dark 
matter mass density profile, while the dot-dashed cyan curve represents a iNavarro. Frenk fc White] l ll996f) mass density profile with 
A/200 = 1.3 x 10 12 Mq, which corresponds to the median halo mass of our sample of simulated galaxies. Right: Median azimuthally- 
averaged surface brightness profile in the V-band for the randomly-oriented simulated galaxies. The solid red curve is the median surface 
brightness profile for all stars, whereas the solid black curve is the median profile for the bulge+halo. The finely (respectively sparsely) 
hatched region encloses the 14 th and 86 th (respectively 5 th and 95 th ) percentiles of the bulge+halo surface brightness profile. The dashed 
green curve represents a power-law profile (functional form given in the legend) with a normalisation chosen to match the median surface 
brightness profile at large radii. No single power-law can match the simulated bulge+halo profiles over all radii. At radii of r < 30 kpc, 
a 'bump' is present, signaling the transition from accretion-dominated to in situ-dominated stars. 



a power-law to the density profile due to accreted stars 
at radii of r ^ 30 kpc, we still find a st eep logarithmic slop e 
of « —2.9 (i.e., similar to that found bv lCooper et al.ll201Cf ). 
The fact that our total (accreted+in situ) profiles show ev- 
idence for an increa se in normalisation at small radii (at 
r ~ 30 kpc), whereas ICooper et al.l (|2010l ) find, if anything, 
a slight flattening, suggests that in situ star formation also 
plays an important role at small radii. We demonstrate that 
this is indeed the case in Section [4] 

It is interesting to compare the mass density /surface 
brightness profiles of the simulated galaxies with those in- 
ferred for local galaxies. We first consider M31, which may 
represent the best test of the models, as its relatively close 
location and our external view point allow observations of 
low surface brightness levels down to ~ 32 mag arcsec -2 and 
hence enable the tracing of the radial profiles out to large 
projected radii. 

In Fig. [4] we compare the V-band surface brightness 
profiles of the simulated galaxies to that derived for M31. 
To date the surface brightness profile of M31 has mainly 
be derived along (or close to) the minor axis. For the pur- 
poses of comparison we therefore also measure the minor 
axis profile of our simulated galaxies. This is done by ori- 
enting the simulated galaxies to an edge-on configuration 
(with x — y plane being parallel to the plane of the disc) and 



measuring the surface brightness as a function of the dis- 
tance z f rom the disc midplane within a cylindrical distance 
(R — y/ x 2 + y 2 ) of 5 kpc from the system center. Unfor- 
tunately, the mass resolution of our simulations prevents us 
from measuring the minor axis profile on a system-by-system 
basif0, so we have stacked all 412 edge-on galaxies to derive 
a mean minor-axis surface brightness profile. This is repre- 
sented by the solid black curve in Fig. [4] For comparison, 
w e show the minor axis surface brightness measurem ent s 
of Pritchet fc van den Berghl l| 19941) . llbata et all (|2007l ). and 
iGilbert et ail (l2009al ). 

Quite remarkably, the shape of the stacked simulated 
profile matches that of the observed surface brightness pro- 
file of M31 over almost two decades in radius (from r ~ 2 
kpc to beyond 100 kpc), including the change in slope that 
occurs at radii of 20 — 30 kpc. 

Note that it is the shape of the profile, rather than its 
normalisation, that is the more critical test of the mod- 
els, since the normalisation of the observed profile requires 
making an uncertain conversion between the tracer pop- 
ulation (Red Giant star counts in this case) to the total 



6 This was possible for the azimuthally-averaged profiles in Fig. [3] 
because of the much larger surface area used for measuring the 
profiles. 
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Figure 4. Comparison of the minor axis V-band surface bright- 
ness profile for the simulated galaxies with that of M31. The solid 
black curve represents the stacked profile derived from all 412 sim- 
ulated disc galaxies. The dotted black curve is the stacked sim- 
ulated profile shifted down by 0.5 mags arcsec -2 . The solid red 
curve represents the de Vaucouleurs profile fit to the minor axis 
profil e of the M31 "bulge" derived by ijPritchet fc van den Berghl 
1994, PvdB). The dashed gre en line represents the power-law fit 
to the outer halo of M31 bv llbata et al.l d2007ll. The solid blue 
circles represent the measurements of lGilbert et al.l l l2009al ). The 
shape of the stacked profile is in very good agreement with that 
derived for M31 over a wide range of radii. 



surface brightness. Typically this is done by assuming that 
the spheroid has an identical stellar luminosity function to 
that of local dwarfs or globular clusters and therefore one 
can use the relationship between integrated light and the 
same tracers in these local systems to re-normalise the halo. 
Of course this procedure will only be accurate if the stel- 
lar halo indeed has a similar stellar luminosity function to 
the local d warfs or globular clusters. For th e data plotted 
m Fig. H IPritchet fc van den Berghl j 19941 ) used the ob- 
servations of the Mil ky Way globular clust er 47 Tuc for 
the con version, while llbata et al.l (|2007l ) and iGilbert et alj 
(2009a) shifted their data vertica lly to map onto the profile 
IPritchet fc van den Berghl (| 19941 ) at radii of ^ 10 - 30 kpc. 
In spite of the potential for important systematics in nor- 
malisation of the observed profile, the stacked profile agrees 
with the data to an accuracy of « 0.5 mags, arcsec -2 . This 
is actually better than one might have naively expected, es- 
pecially given that the system-to-system scatter in the sim- 
ulated profiles is w 3 mags, arcsec -2 (see Fig. [3]). 

We turn now to the Milky Way. Bearing in mind 
that observations of the extended stellar distribution of 
the Milky Way are typically limited to r < 40 kpc, 
a wide variety of studies have found logarithmic slopes 
for the density profile of the stellar halo ranging from 
from —2.5 to —3.75, dep ending on the data and m e thod- 
ology that were used ilHarrlsl Il976|; IZinnl Il985l; ISaha 



1996 


; 


Ivezic et al. 200d; Morrison et al. 


l200d; lYannv et al. 


200C 


Chiba & Beers 20001: Chen et al 


200 ll; Siegel et al. 


2002 


;l Vivas & Zinnll2006; Bell et al. 


200S 


June et al.ll2008l; 


Smith et alj 120091; ISesar et al. 2011 


; Deason et al.l l2011bl). 



ology that were used lliiarrisl 119/rj; IZiinnl 119851; ISaha 
Il985l : iPreston. Shectman Beersl Il99ll ; IWetterer fc McGraw 



This is quite compatible with results of the cosmological 
hydrodynamical simulations. Interestingly, some of the lat- 
est observations of the Milky Way halo, which are based 
on large samples of stars of the SDSS, find evidence for 
a change in the slope densi t y profile at radii of around 
» 30 kpc (ICarollo et al.ll2007l ; fjuric et al.ll2008l ; ISesar et all 
I2011I ; Id eason et al.ll2011h] r which is remarkably consistent 
with our simulations. The same studies also find strong evi- 
dence that the inner halo has an oblate distribution (rather 
than spherical), with a minor/major axis ratio of 0.5 — 0.7. 
We present a study of the shape and kinematics of the simu- 
lated galaxy haloes in McCarthy et al. (in prep.). However, a 
visual inspection of the simulated azimuthally-averaged and 
minor-axis profiles in Fig. [3] and 3] already reveals that our 
spheroids are indeed flattened within 30 kpc or so. 

Observations of stellar haloes of other (more distant) 
external galaxies are more difficult. Currently, only a few 
stellar haloes of individual galaxies have been detected, and 
only down to levels of ~ 28—29 mag arcsec -2 (|de Jong et alj 
2007), corresponding to the inner regions. One can go deeper 
by st acking distant galaxies. For example, IZibetti et al.l 
(2004) stacked more than 1000 stellar haloes of edge-on spi- 
ral galaxies observed in the SDSS and found that, out to 
about 25 kpc (Ey- ~ 31 mags, arcsec -2 ), the average den- 
sity profile can be fitted relatively well with a power-law 
p, ~ r -3 , very similar to what the simulations predict. 

With our large sample of simulated galaxies (~ 400), 
we are in an excellent position to characterise the system- 
to-system scatter in the normalisation of the p* and £y- 
profiles, which is another potential test of the simulations. 
The solid black curves represent the median profiles of the 
bulge+halo components only. The finely-shaded regions en- 
close the 14 th and 86 th percentiles (i.e., la) and the sparsely- 
shaded regions enclose the 5 th and 95 th percentiles (2cr) 
of the stellar bulge+halo profiles. The 2er scatter in p„ is 
approximately an order of magnitude over the range of 
10 — lOOkpc, corresponding to a scatter in £y of ~ 3 
mags. /arcsec 2 over the same range. [In the Appendix, we 
test the numerical convergence of the stellar mass density 
profiles. We demonstrate there that our results (both median 
and scatter) are quite robust to a factor of eight increase in 
the mass resolution.] 

It is also worth commenting briefly on the structure 
of the dark matter halo. In addition to the lines corre- 
sponding to the stellar components described above, the 
left panel of Fig. [3] also shows the median dark matter den- 
sity distribution (as the blue cu rve) and an analytic NFW 
l|Navarro. Frenk fc White] Il996h profile (the cyan curve). 
The analytic halo has a mass M200 = 1.4 x 1O 12 M0 (which 
corresponds to the median halo mass of our sample of simu- 
lated galaxies) and a concentration parameter of C200 = 9.4, 
which places it on the M200 — C200 relation measured by 
iGao et alj l|2008h for the Millennium Simulation. The dot- 
dashed cyan curve is not a fit to the dot-dashed blue curve. 

A comparison of the dot-dashed blue curve with the 
solid red curve demonstrates that dark matter dominates 
over the stellar contribution down to very small radii 
(the two become comparable at a few kpc). Comparison 
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Figure 5. Left: Median spherically-averaged stellar metallicity profiles. The solid red curve is the median stellar metallicity profile for 
all stars, whereas the solid black curve is the median stellar metallicity profile for just the bulgc+halo component. The finely (sparsely) 
hatched region encloses the 14 th and 86 th (5 th and 95 th ) per centiles of the stellar halo metallicit y profi l e. The solid blue, green, cyan, 
and orange circles represent M31 metallicity measurements of lGilbert et al.l j2009al ), iKalirai et al.1 (2006), Richar dson et al] {2009), and 
iKoch et all l|2008r i. respectively. The observational data have been scaled to the same set of solar abundances. Right: Median metallicity 
distribution function (MDF) in three radial bins. Significant negative metallicity gradients are a ubiquitous feature of the simulated 
galaxy populations. 



of the dot-dashed blue and cyan curves shows that the 
NFW analytic form provides an excellent description of 
the dark matter density profile beyond about r ~ 10 
kpc (i.e., ~ 0.05r2oo)- Within this radius, the dark mat- 
ter profile steepens, presumably as a result of contrac- 
tion of the d ark matter halo due to the condensation of 
baryons (e.g.. iBlumenthal et al.l Il98rj ; iGnedin et ail 12004 
iDuffvet al.ll2010r i. 



Given that the NFW profile drops off as ~ r~ 2 at in- 
termediate radii and as ~ r~ 3 at large radii, this implies 
that the stellar mass density profile drops off significantly 
faster than that of the dark matter, as is clearly visible in 
the left panel of Fig. [3] This is also consistent with the 
findi ngs of previous theoretical studies of the st ellar halo 
(e.g. iBullock fc Johnstonll2005l ; lAbadi et al.ll2006l ). The dif- 
ference in the stellar and dark matter distributions likely 
stems from the fact that the stars in infalling satellites are 
more tightly bound, and thus less susceptible to tidal disrup- 
tion, than their dark matter haloes. In addition, the most 
massive satellites (which also have the highest stellar mass 
fractions) will spiral more quickly into the center due to 
dynamical friction, which scales as mass squared. Finally, 
in situ star formation will occur preferentially in the cen- 
tral regions (due to the higher gas densities there). These 
effects would produce a more centrally-concentrated stellar 
distribution, as seen in the simulations as well as in the ob- 
servations. 



3.2.2 [Fe/H] and [a/Fe] radial profiles 

The radial distribution of stellar metallicity is a particu- 
larly interesting test of the simulations, since previous sim- 
ulations/models in the ACDM context have not been able 
to reproduce the significant gradients seen in M31 and the 
Milky Way. It is not presently clear why this is the case. It 
may signal a fundamental problem with galaxy formation in 
the ACDM context or that previous models of stellar haloes 
neglected important physics, such as in situ star formation 
or a sufficiently realistic implementation of chemodynamics. 

In the left panel of Fig. [5] we show the me- 
dian spherically-averaged metallicity profile and system-to- 
system scatter about the median. The metallicity, [Fe/H], 
of each star particltQ is computed by taking the logarithm 
of the ratio of the iron-to-hydrogen mass fractions of the 
star particle and subtracting the logarithm of the iron-to- 
hydrogen mass fraction ratio of the Sun (which we assume 
is 0.00156, co nsistent with the rece nt measurements and 3D 
modelling of lAsplund et all [2005). A spherically-averaged 
metallicity profile is computed for each galaxy by comput- 
ing the mean of metallicity (i.e., <[Fe/H]>) of star particles 



7 For consistency with the way the chemodynamics is done in the 
simulations, we used so-cal led 'smoothed m etallicitics' rather than 
'particle metallicitics'. See Wiersma et al. (2009b) for discussion. 
Switching to 'particle metallicities' has the effect of shifting the 
profiles in the left panel of Fig. [5] down by 0.2 dex, which, 
as discussed later in the paper, is smaller than the uncertainties 
in the adopted empirical nucleosynthesis yields and Type la SNe 
rates. 
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in radial bins. The median profiles plotted in Fig.[5]are then 
computed by taking the median of all of individual metal- 
licity profiles in radial bins. 

Interestingly, the stellar metallicity distribution shows a 
prominent negative gradient over all radii, with the steepest 
decline exhibited over the range 20 kpc < r < 40 kpc. This 
is independent of whether the disc component is included 
or not. The solid black curve and shaded regions indicate 
that negative metallicity gradients are a ubiquitous feature of 
the simulated galaxy populations. (This statement is further 
confirmed in the middle panel Fig. [9] in Section [4j which 
shows system-to-system scatter in the difference in the mean 
metallicity within and beyond 30 kpc.) 

The solid blue, green, cyan, and orange cir- 
cles in Fi g. [5] represen t metal lic ity measurem e nts o f 
M31 from IGilbert et all (l2009all. IKalirai et al l l|2006h . 
IRichardson et all l|2009ll and iKoch et all |2008f) . respec- 
tive!^. The simulated profiles have a rem a rkably simi- 
lar shape to that derived by iGilbert et all (|2009al ) and 
IKalirai et all (2006), which shows a drop of w 0.7 — 0.8 dex 
over the wide radial range of ~ 10 — 150 kpc. The sim- 
ul ated profiles are also co nsistent with the measurement s 
of IRichardson et all l|2009Il (see also lChapman et af1 [2006l 
which pr obe a much s malle r dynamic range. The profile de- 
rived by Koch et al . (2008), however, is somewhat steeper 
than the predicted prof iles, particula r ly at v e ry large radii. 
It is w orth noting that IGilbert et all (|2009al ). IKalirai et all 
l|2006h . and IRichardson et all (|2009h all employed similar 
methods in deriving the metallicity (in particular, they used 
spectroscop i cally- calibrated photometric methods), whereas 
IKoch et al. (2008) used a new spectroscopic method that de- 
rives abundances from the calcium triplet (CaT) part of the 
spectrum. It is not presently clear what the origin is of the 
difference between the different data sets. 

It is also encouraging that the normalisation of the sim- 
ulation metallicity profiles agrees as well as it does (to within 
0.2-0.3 dex accuracy) with the observations of M31. The sim- 
ulations adopt empirical nucleosynthesis yields and Type la 
SNe rates, both of which a re uncertain by a factor of a few 
(see IWiersma et all l2009bh , implying that there is freedom 
to vertically shift the simulated profiles up or down by a few 
tenths of a dex. Another potentially relevant factor is the se- 
lection of the observational fields. Often t hey are sele c ted to 
have a relatively high sur face brightness. iFont et all (120081 ) 
and IGilbert et all l|2009al ) have shown that surface bright- 
ness tends to be positively correlated with metallicity, so it 
may well be the case that the current measurements are bi- 
ased somewhat high. Large panoramic observations of M31 
(such as PAndAS) will soon shed light on whether this bias 
is important. 

For the Milky Way, metallicity measurements out to 
large galactocentric radii are much more challenging. How- 
ever, recent results based on SDSS observations suggest 
suggest a similar drop of ~ 0.7 — 0.8 dex from the Solar 



8 Where appropriate we have adjusted the observational metal- 
licity measurements (by applying a constant offset) to account 
fo r differences in the a ss umed solar abundanc es in the studies 
of IGilbert et alj d2009all IKalirai et aP l|2006l l. IRichardson et al] 
d2009lV and | Kochet*ljT i2008f ) from the solar abundances of 
lAsplund et al .1 l|2005l l. which we use to normalise the simulated 
profiles. 



neighbourhood out to « 30 — 40 kpc (ICarollo et al. 
Ilvezic et all 120081: ICarollo et all I2OI0I ; Ide Jong et all 
but see ISesar et al.ll201ll l. 



2007 



2010, 



As alluded to above, most previous theoretical stud- 
particular ly those based o n accretion-only hybrid 



models (e.g IFont et al l l2006al ; |Pe Lucia fc Helmil 120081 ; 



ICooper et akl feoiO). did not find significant gradients. To 
our knowledge, previous studies based on hydrodynamic 
re-simulations of disc galaxies (e.g., lAbadi et all 120061 ; 
IZolotov et alll2010h did not examine the radial distribution 
of the stellar metallicitjjf]. As we demonstrate in Section [4] 
the change in slope at r ~ 30 kpc marks the transition from 
accretion (r > 30 kpc) to in situ (r < 30 kpc) stellar halo 
formation. 

While the spherically-averaged profiles of the simulated 
galaxies display prominent gradients, at any given radius 
there is a wide range of stellar metallicities. To illustrate the 
spread in metallicities as a function of radius, we show the 
median stellar metallicity distribution function (MDF) for 
three radial bins in the right panel of Fig. [5] For all three bins 
the MDF is similar in shape (in particular, fitting a Gaus- 
sian distribution yields similar widths of a ~ 0.5 dex), but 
the peak shifts to lower metallicities with increasing radius, 
which is what gives rise (mathematically at least) to the 
metallicity gradients in the left panel of Fig. [5] A Gaussian 
provides a reasonable match to the MDFs near the peak and 
for higher metallicities, but it systematically under-predicts 
the fraction of halo stars with low metallicities. The MDFs 
in the right panel of Fig. [5] are qualitatively similar to that 
found for the M31 ste llar halo (Gilbert et al . in prep) and 
other nearby galaxies (|Mouhcine et al.ll2005bT) . 

Since the simulations include enrichment from both 
Type II and Type la supernovae, as well as stellar mass 
loss (i.e., AGB stars), we can also make predictions for the 
[a/Fe] radial profiles, which should be measurable in galaxies 
like M31 in the near future. The so-called a elements (such 
as oxygen, magnesium, and silicon) provide an additional 
constraint, as they are produced primarily in massive stars. 
In the left panel of Fig. [6] we show the median spherically- 
averaged [a/Fe] profile and scatter, where we have used both 
oxygen and magnesium to represent a. We have normalised 
the profiles using the oxygen-to-iron and magnesium-to-iron 
mass fractions of th e Sun (5.2887 and 0.5355, respectively, 
lAsplund et al.ll2005f ). Interestingly, while the metallicity (i.e. 
[Fe/H]) profile shows a strong radial dependence (Fig. [S}, 
the [a/Fe] distribution does not: [a/Fe] rises by only ~ 0.1 
dex from the center to the outskirts. This implies that the 
physical mechanism responsible for producing the iron does 
not change appreciably with galactocentric distance. 

It is interesting to elucidate the lack of a strong ra- 
dial dependence in the [a/Fe] profiles presented in the left 
panel of Fig. [5] The GIMIC simulations track not only the 
total iron mass fraction, but also the iron mass fraction pro- 
duced by Type la supernovae only. This therefore allows us 
to determine, for any given star particle, the individual con- 
tribution of Type la to the iron mass fraction. In the right 



9 But note that hydrodynamic re-simulations have been used to 
study the radial distribution of metallicity in ea rly-type galaxies 
(e.g.. lKobavashill2004 |2005|; iGibson et al.ll2007l) and significant 
gradients were shown to be common in these systems as well. 



Stellar haloes of disc galaxies 13 




Figure 6. Left: Median spherically-averaged stellar [ct/Fe] profiles. The solid black (light blue) curve is the median [O/Fe] ([Mg/Fe]) 
profile for the bulge+halo component. The finely (sparsely) hatched region encloses the 14 th and 86 th (5 th and 95 th ) percentiles. The 
dashed black (light blue) curve shows the effect on the median [O/Fe] ([Mg/Fe]) profiles if the rate of Type la SNe were to be increased 
by a factor of 3 in the simulations. Right: Median spherically-averaged stellar [Fe/H] profiles for the spheroidal component using the total 
mass of iron (green) and just the iron produced by Type la supernovae (orange). The finely (sparsely) hatched regions enclose the 14 th 
and 86 th (5 th and 95 th ) percentiles. There is no significant radial variation in the simulated [a/Fe] profiles, owing to the dominance of 
metal production by Type II SNe at all radii. 



panel of Fig. [S]we show the [Fe/H] profiles for the total iron 
mass fraction and for the iron produced solely by Type la 
supernovae. As it can be clearly seen, Type la supernovae 
make only a minor contribution to the metallicity of the 
stellar halo. Type II supernovae are therefore the dominant 
producers of metals at all radii. The very mild gradient in 
[a/Fe] is due to the relatively larger importance (though still 
small in an absolute sense) of Type la for stars located in 
the inner regions. 



It is important to acknowledge potential caveats in the 
above prediction. Besides the nucleosynthesis yields, the rate 
of Type la supernovae is very uncertain. Rather than try 
to calculate the rate of Type la supernovae using a phys- 
ical model, the simulations adopt empirical Type la rates, 
which are still uncerta in by a factor of a few (see Fig. A6 of 
IWiersma et al.|[2009bh . If the Type la rate is actually a fac- 
tor of a few higher than current estimates imply, this would 
have the effect of producing a steeper positive [a/Fe] gradi- 
ent (see the dashed curves in the left panel of Fig. [5] which 
show the effect of boosting the number of Type la SNe by 
a factor of 3), since the stars in the central regions were 
formed more recently (see Section Q and therefore the gas 
from which they formed had more time to be polluted by 
Type la supernovae. 



4 PHYSICAL ORIGIN OF THE RADIAL 

DENSITY AND METALLICITY PROFILES 
OF THE STELLAR HALOES 

In the previous section we showed that the GIMIC simulations 
produce disc galaxies with stellar haloes that bear a strong 
resemblance to those observed around nearby galaxies, in- 
cluding M31 and our own Milky Way. In the present section 
we use the simulations to trace back the physical origin of 
the stellar haloes. We first quantify the relative importance 
of accretion vs. in situ star formation for the formation of 
the stellar halo. We show that the break in the mass density 
and surface brightness profiles near r ~ 30 kpc and the sharp 
rise in the stellar halo metallicity profile at similar radii are 
all caused by the transition from accretion-dominated to in 
situ-dominated star formation (Section [4.2p . We also explore 
the time of star formation and assembly of the stellar halo 
(Section 14. 3p and the degree of radial migration of the in 
situ component (Section 14.41) . 



4.1 Distinguishing accretion and in situ formation 

To assess the relative importance of accretion vs. in situ 
star formation we must first define what we mean by these 
terms. We define in situ star formation as star formation 
that takes place in the most massive subhalo of the most 
massive progenitor (MMP) FoF group of the present-day 
system. Stars are said to be 'accreted' if they formed either 
in another FoF group or in a satellite galaxy (i.e., not the 
dominant subhalo) of the MMP. In practice, only a small 
proportion (2%) of the stellar halo is from stars that formed 
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in satellites of the MMP and were then stripped; most of the 
accreted stars were formed in other FoF groups (i.e., other 
haloes) before they became satellites of the MMP. 

The procedure by which we identify the MMP and de- 
termine whether a given star particle was accreted or formed 
in situ is as follows. For a given Milky Way-mass system at 
z = (lwe select all of the dark matter particles within r2oo 
belonging to the dominant subhalo. We use the unique IDs 
assigned to these particles to identify the FoF group in pre- 
vious snapshots that contains the largest number of these 
particles. This FoF group is said to be the MMP and the 
dominant subhalo of that FoF group is assumed to be the 
most massive progenitor of the dominant subhalo of the sys- 
tem at z = 0. Each of the star particles that comprise the 
z = stellar halo is tracked back in time to the FoF group 
and subhalo it belonged to, if any, when it formed from a gas 
particle. If, at the time of star formation, the star particle 
is a member of the dominant subhalo of the MMP, the star 
particle is said to have formed in situ. If it formed in another 
FoF group or in a non-dominant subhalo of the MMP, then 
it is said to have been accretecP^I. 

The number of snapshots output for each of the GIMIC 
volumes is approximately 60, and they span the redshift 
range z = 20 to z = 0. The time of star formation of each 
star particle will, in general, not correspond to the redshift 
at which the snapshot data was written. We therefore must 
use the snapshot that is closest in time to identify where 
the star particle was when it formed. This could, in certain 
circumstances, lead to a misidentification of the location of 
star formation (i.e., to which subhalo it belonged when it 
formed). However, we have explicitly checked that this is not 
an important issue by re-running our tracing code using only 
every second snapshot. We find that the in situ mass fraction 
of our simulated galaxies calculated using only every second 
snapshot agrees with that inferred using all the snapshots 
to an accuracy of 1.6% on average. We therefore conclude 
that the redshift sampling of our snapshot data is sufficient 
to deduce whether a given particle formed in situ or was 
accreted. 

We trace each of the z — stellar halo particles back 
in time as far as z = 10 to determine whether they formed 
in situ or were accreted. We ignore the very small fraction 
(0.01%) of halo stars that formed before z = 10. 



4.2 Stellar mass density and metallicity profiles 
by formation mode 

In Fig. [7]we show the median spherically-averaged mass den- 
sity profile of systems split according to formation mode 
(i.e., in situ vs. accreted). A comparison of the dot-dashed 
red and dashed blue curves demonstrates that in situ star 
formation contributes significantly to (dominates) the stel- 
lar halo within a radius of 30 (20) kpc. Beyond 30 kpc the 
majority of the stars were formed in satellites prior to accre- 
tion. This transition accounts for most of the deviation of 



10 We have also allowed for the possibility that the star particle 
may have formed outside of bound substructures but find that 
this mode of formation is rare, accounting for 2% of the stellar 
halo by mass. 
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Figure 7. Median spherically-averaged stellar mass density pro- 
files. The dotted and solid black curves represent all stars and 
stars in the spheroid only, respectively. The dotted black curve 
represents all stars (disc+spheroid). The dashed blue curve rep- 
resents accreted halo stars (i.e., stars that formed outside the 
most massive progenitor). The dot-dashed red curve represents 
spheroid stars that formed in situ. In situ star formation con- 
tributes significantly to (dominates) the stellar halo within a ra- 
dius of 30 (20) kpc. 
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Figure 8. Median spherically-averaged metallicity profiles. The 
curves have the same meaning as in Fig. [7] For r < 50 kpc, stars 
formed in situ have higher metallicity and there dominance of the 
stellar halo by mass at small radii is what produces the strong 
radial variation in [Fe/H], 
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Figure 9. The current stellar mass fraction of the bulge+halo component due to in situ star formation as a function of (Left) mean 
metallicity, [Fe/H], (Middle) the metallicity difference between the inner and outer bulge+halo, A [Fe/H] (calculated within and beyond 
30 kpc), and (Right) system formation epoch, 2f orm (i.e., the redshift when half the present-day total mass (gas+stars+DM) is in place). 
Each data point represents a galaxy. Galaxies with higher in situ mass fractions generally have higher mean metallicities, more prominent 
metallicity gradients, and are dynamically older, although the system-to-system scatter in these trends is large. 



the total stellar halo mass density profile from a pure power- 
law distribution of r~ 3,5 at small radii (compare the solid 
black and dashed green curves in the left panel of Fig. [31). 

Averaged over all of the simulated galaxies, in situ stars 
account for approximately 68% of the current mass of the 
stellar spheroids (with considerable system-to-system scat- 
ter, as shown in Fig. 0), while accreted stars constitute ap- 
proximately 32%, with 28% of the mass formed in satel- 
lites prior to accretion, 2% formed in satellites after ac- 
cretion, and 2% accreted smoothly. Our inferred average in 
situ fraction i s there fore larger than that found recently by 
IZolotov et all l|2009l ), who found that 3 out of their 4 re- 
simulated galaxies had in situ mass fractions of < 20%. How- 
ever, part of this difference is likely due to the fact that we 
have made no attempt to distinguish between bulge and halo 
components (see Section [2.11 for a discussion of this p oint). 
Consistent with this is that both lAbadi et alj |2006l ) and 
lOser et "all |2010T ) find in situ mass fractions of ~ 40 — 50% 
and, like us, these authors have not decomposed the bulge 
from the hak£E 

In terms of our simulations, if we limit our 
comparison to r > 20 kpc (well beyond any likely bulge com- 
ponent), the in situ fraction drops to 20% while the accreted 
fraction rises to 80% (with 68% formed in satellites prior to 
accretion, 3% formed in satellites after accretion, and 9% 
accreted smoothly). 

In Fig. [8] we show the median spherically-averaged stel- 
lar metallicity profile split by formation mode. Here we 
see that the gradient in the total stellar metallicity profile 
(solid black curve) is driven by the generally high metal- 
licity of stars formed in situ combined with the fact that 
the in situ stars dominate the stellar halo by mass in the 
inner regions. The metallicity of accreted stars shows no 



11 On e possible reason why lAbadi et al. (2006) and O ser et al. 
i20irj ) find slightly lower in situ mass fractions than us, is that 
their simulations ignore metal-line cooling, which is expected to 
be important. 



strong variation with radius (except possibly in the cen- 
tral few kpc). The fact that the overall metallicity gradi- 
ent is established as the result of a transition from in situ- 
dominated to accretion-dominated stars, appears to resolve 
the longstanding pro blem of the inability of purely accretion- 
based models (e.g . . | Font et al.l l2006ah |Pe Lucia fc Heliml 
I200SI ; ICooper et alj|2010h to reproduce the strong variation 
in th e mean stellar meta l licity with radius observed in M31 
(e.g., iKalirai et all 120061; iKoch et all l2008h and t he Milky 
Way (e.g- ICarollo et alj|2007tlde Jong et alj|2010l ). 

In Fig. [9] we show scatter plots of the in situ mass frac- 
tion as a function of the mean metallicity of the bulge+halo 
(left panel) and the metallicity difference between the mean 
metallicity within and beyond 30 kpc (middle panel). Each 
data point represents a simulated galaxy. Although there 
is considerable system-to-system scatter, there are identifi- 
able correlations between the in situ mass fraction and the 
mean metallicity and metallicity difference, in the sense that 
galaxies with higher in situ mass fractions have higher mean 
metallicities and more prominent metallicity gradients. We 
discuss the right panel of Fig. [9] below. 



4.3 The time of star formation and assembly of 
the stellar halo 

We now examine the age of the stars and the assembly time 
of the spheroidal stellar component. In Fig. \W\ we show 
profiles of the median spherically-averaged formation look- 
back time (i.e., the age) of stars, split by formation mode. 
Accreted stars (dashed blue curve) are typically quite old 
(~ 11 — 12 Gyr), whereas most of the stars that formed 
in situ and which are located within 30 kpc or so, have a 
median age of w 7 — 8 Gyr. 

The finding that stars that were formed in situ are, 
on average, younger than accreted stars, is consistent with 
our expectations based on the hierarchical growth of struc- 
ture. In particular, stars that formed in situ are, by defini- 
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Figure 10. Median spherically-averaged stellar age (rf orm ) 
profiles. The solid black curve represents all stars in the 
spheroidal component. The dotted black curve represents all 
stars (disc+spheroid). The dashed blue curve represents accreted 
stars in the spheroid. The triple dot-dashed red curve represents 
spheroid stars that formed in situ. The dot-dashed blue curve 
represents the time of accretion of accreted stars (i.e., the time 
when the system in which the stars formed became a satellite of 
the MMP). The accreted component is old (fa 11 - 12 Gyr), but 
the lookback time at which it was assembled varies strongly with 
radius, with the outer most regions having been accreted quite re- 
cently. The in situ component is younger and was typically formed 
at z S3 1 - 1.5. 



tion, those that formed in the MMP, whereas accreted stars 
formed in less massive and typically dynamically older satel- 
lites. For a 1 - 4 x 10 12 M Q (DM) halo the epoch of forma- 
tion, defined as the r edshift where half th e final dark mat- 
ter mass is in place dLacev fc Cold fl993T l . is z « 1 — 1.5 
(e.g., IWechsler et al.ll2002T l. which corresponds roughly to 
the average star formation time of the stars formed in situ. 
Of course, baryons feel not only the effects of gravity, but 
also those of non-gravitational processes such as cooling and 
heating due to feedback (e.g., from supernovae). Feedback 
in particular can result in the delay (or prevention) of signif- 
icant star formation. C09 have examined the star formation 
histories of galaxies in the GIMIC simulations and found that 
the peak occurs at z ~ 1 — 3 for normal galaxies, but at 
z ~ 3 — 6 for dwarf galaxies (which are the systems in which 
the vast majority of the accreted stars formed). 

Returning to Fig. [9] we plot in the right hand panel 
the in situ mass fraction of the galaxy against the formation 
epoch of the system [i.e., the redshift for which the MMP had 
a total mass due to gas+stars+DM that is half of the final 
[z — 0) total mass]. There exists a correlation (again with 
considerable scatter) between in situ mass fraction and for- 
mation epoch, in that older systems have higher in situ mass 
fractions. Physically this makes sense, as older systems have 
been able to form stars quiescently for a longer period of 
time, while dynamically younger systems accreted relatively 



larger amounts of mass recently. The fact that dynamically 
older systems have high er in situ mass frac tions is consis- 
tent with the findings o f lZolotov et ail < |2009h . although our 
improved statistics shows that the scatter in this correlation 
is large. 

It is interesting that, even though the accreted stars 
are quite old (as indicated by the dashed blue curve in 
Fig. 1101) . many of them were accreted recently, as indicated 
by the dot-dashed blue curve in Fig. [TUJ In agreement with 
hybri d simulations of pure acc r etion-based halo f ormation 
(e.g., iBullock fc Johnston! liooBI ; iFont et all [200681 ). we find 
that the inner accreted halo was formed a long time ago 
(Tacc ~ 9 — 10 Gyr), but that the outer accreted halo was 
assembled as recently as 5 Gyr ago. We therefore expect to 
see increased spatial lumpiness and more coherent velocity 
structure at large radii, since there has been less time since 
accretion for phase mixing and violent relaxation to act. We 
plan to examine the degree of substructure in the spheroidal 
component in a future paper. 



4.4 Radial migration of in situ stars 

We have shown that it is in situ star formation that is re- 
sponsible for the change in the nature of the stellar halo at 
r < 30 kpc. Did these stars for m at small radi i and w here 
the n dispersed, as sug gested by IZolotov et alj (|2009h (see 
also 



Purccll et al. 2010)? 
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Zolotov et al. (2009) found that the in situ stars in their 



re-simulations were formed from gas that was brought into 
the main system in smooth 'cold flows'. In their simulations 
the stars formed at very small galactocentric distances of 
~ 1 kpc. These stars were then dispersed out to larger radii 
(as far as 20 kpc or so) by z ~ 2, with major mergers acting 
as the primary dispersal mechanism. 

In Fig. [TT] we show the distribution of initial (i.e., 
at the time of star formation, Zf orm ) and final (z — 0) 
radii of stars that formed in situ in the GIMIC simulations. 
This distribution is computed by first creating a regular 
(50 x 50) 2D grid in initial and final radii (in log space, 
from log 10 (r) = 0.0 — 2.5 in steps of 0.05) and then count- 
ing the number of in situ star particles in each pixel. The 
map is normalised by the total number of in situ star par- 
ticles. The red, yellow, green, and blue iso-density contours 
enclose 50%, 68%, 90%, and 95% of the total number of in 
situ star particles, respectively. Fig. [TT] represents a stacked 
distribution of all the in situ stars of all simulated galaxies. 

Fig. [TTJ shows that the vast majority of in situ stars 
are presently located at, and were also formed at, relatively 
small radii (r < 30 kpc). At these radial distances there is 
a relatively large scatter about the r[z = 0] = reform] re- 
lation, indicating that in situ stars have migrated in both 
directions (i.e., toward and away from the present galaxy 
center). A very small fraction of the in situ stars were formed 
at relatively large radii (r > 30 kpc). According to Fig. [TTJ 
these stars have since migrated inward somewhat. Naively, 
it seems surprising that star formation can occur at such 
large galactocentric radii. The density and temperature cri- 
teria imposed in the simulations for placing gas particles on 
the star-forming EOS (see Section 2) guarantee that these 
in situ stars formed from gas that was both dense and cold. 
We have analysed the in situ star formation in the outer 
halo in depth. First, such star formation is virtually absent 
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Figure 11. Distribution of initial (i.e., at redshift of star forma- 
tion, 2f orm ) and final (z = 0) radii of stars formed in situ. The red, 
yellow, green, and blue filled iso-density contours enclose 50%, 
68%, 90%, and 95% of the total number of in situ star particles, 
respectively. The dotted curve represents the r[z = 0] = r [zf orm ] 
relation. On average, in situ stars formed at radii comparable 
to where they are presently located, although migration in both 
directions (towards and away from the center) is significant. 



in the simulated galaxies at the present day. Like the ma- 
jority of the population, these stars formed at z > 1. This 
suggests a physical (rather than spurious numerical) origin 
to these stars. Second, stars formed beyond r[zf orm ] > 20 
kpc show up as 'streaks' in a plot of r[zf orm ] against Zf or m 
(i.e., collections of stars formed at very similar time over a 
wide range of 3D radii), whereas at small radii there is a 
larger and more homogeneous spread in ages. This suggests 
that the large-radii stars are forming in streams (filaments) 
of dense/cold gas. Finally, examination of simple movies of 
the high-resolution simulation shows that the star forma- 
tion indeed occurs in filaments of cold gas, which has re- 
cently been ram pressure-stripped from infalling satellites. 
This type of star for mation has been seen i n SPH simula- 
tions previously (e.g.. [Puchwein et al.|[201ol ). For idealised 
test cases it has been shown that the SPH for malism can in- 
hibit the development of fl uid instabilities (e.g. jAgertz et all 
2007:; iMitchell et alj|2009h . which could potentially disrupt 
the clouds/filaments. Thus, it is possible that our simula- 
tions overestimate the degree of in situ star formation at 
such large radii. In any case, these stars represent only a 
small fraction of the in situ population and do not drive 

trends seen in, e.g., Figs. 3-5. 

Our results differ from those of lZolotov et alj l|2009l ) in 
terms of the radius of in situ star formation, with the in 
situ stars in our simulations forming at radii comparable to 
where they are presently located. Our in situ st e llar p opula- 
tion is also younger than that of lZolotov et al. (2009) (who 
found Zform « 3), having formed primarily at z ~ 1 — 1.5, 
which is similar to that also reported recently bv lOser et al.l 



j20ld ) for galaxies with masses similar to those in our sam- 
ple. The differences in the formation radii a nd ages of the 
in situ stars in our simulations and those of IZolotov et al.l 
(2009) would be expected if the feedback in their simulations 
is relatively less efficient in low-mass haloes. This indeed ap- 
pears to the their simulations significantly overpro- 
duce the number of bright satellite galaxi es orbiting Milky 
Way- mass galaxies today (see Fig. 4 of iGovernato et al] 
120071 ). 

More importantly, the above differences imply a differ- 
ent formation scenario for in situ sta r s in o ur simulations 
compared with those of lZolotov et al.l l|2009h . In McCarthy 
et al. (in prep.), we will show that a large fraction of the in 
situ stars formed from the cooling of hot gas ( rather than 
cold flows, as advocated bv lZolotov et al]|2009l ) during the 
initial assembly of the disc. This has profound implications 
for a wide range of properties of these stars today, including 
their spatial distribution, kinematics and metallicities. 

4.5 A note on "bulgeless" galaxies 

iKormendv et all (|2010l ) have recently demonstrated, using 
high quality photometric and spectropscopic observations, 
that many of the largest nearby disc galaxies lack obvi- 
ous classical bulge components (i.e., are "bulgeless"). These 
authors have argued that it is difficult to reconcile a sig- 
nificant population of bulgeless galaxies with the standard 
hierarchical model for galaxy formation, since mergers are 
a key ingredient of this model and are believed to produce 
bulges. However, recent high resolution cosmological simula- 
tions ha ve been able to successf ully produce bulgeless dwarf 
galaxies |Governato et al1l2010l ) . The key to this success ap- 
pears to be the implementation of efficient supernova feed- 
back, which preferentially removes the lowes t angular mo- 
mentum gas via outflows (|Brook et alj|201ll ). What about 
more massive simulated galaxies? It turns out that the GIMIC 
simulations produce a number of massive "bulgeless" galax- 
ies, which will be examined in detail in Sales et al. (in prep). 

For the present study, it is interesting to ask whether the 
simulated galaxies with small (fractional) bulge components 
have stellar haloes at all and, if so, what is the contribution 
of in situ star formation. It is of interest to note that real 
bulge less dwarf galaxies do indeed have faint, old stellar ha- 
los (|Stinson et aljj 2009) and so d o the simulated bulgeless 
dwarf galaxies of lGovernato et all (|2010T ). 

In the left panel of Fig. [T^] we compare the median 
spherically-averaged stellar mass density profiles of the total 
bulge+halo and the in situ bulge+halo components for our 
entire sample of galaxies with a subset of galaxies with high 
D/T and a subset with low D/T. Interestingly, we find that 
the total and in situ haloes have larger mass densities in 
systems with high D/T ratios. This is the opposite of what 
one might have naively expected. The reason for this trend 
is provided in the right panel of Fig. 1121 which is a scatter 
plot of D/T versus the total stellar mass (disc+bulge+halo) 
for all the simulated galaxies. As can be seen, the highest 
D/T galaxies have larger total stellar masses (and larger 
stellar mass fractions) on average than galaxies with low 
D/T. This explains why the absolute mass densities of the 
in situ bulge+halo component in high-D/T systems is larger 
than that of this component in low-D/T systems. It should 
be noted that this comparison is being made at essentially 
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fixed total system mass (M200 varies by only a factor of 4 
over our sample). By comparison, the total stellar mass M» 
varies by factor of « 30 across the sample. 

Thus, based on results in Fig. [12] "bulgeless" galaxies 
may actually have the most prominent stellar haloes at fixed 
total system mass, which is an interesting prediction of the 
model that should be testable in the future. The origin of 
the relatively strong D/T-A4* trend over the narrow range 
of total system masses considered in this study is unclear 
and requires further examination. Furthermore, whether or 
not this trend is consistent with observations should be in- 
vestigated. However, as discussed in Section 3.1, care must 
be taken to measure D/T in the same way for both real and 
simulated galaxies. 



5 SUMMARY AND DISCUSSION 

We set out to address the question of the possible dual- 
nature of the extended stellar distributions around the Milky 
Way, M31, and other ~ L* nearby galaxies. In particular, 
what are the relative contributions (by mass) to the stellar 
spheroid from tidal stripping/disruption of satellites and 'in 
situ' star formation? And what is the origin of the metallicity 
gradients observed in the Milky Way and M31? 

We have used the GIMIC suite of cosmological hydrody- 
namical simulations to address these important questions. 
GIMIC is ideal for this purpose as it contains large numbers 
(~ 400) of Milky Way-mass galaxies simulated at relatively 
high- resolution (M g ~ 1.16 x 10 7 /i _1 Mq and a gravitational 
softening of l/i -1 kpc) and with sophisticated prescriptions 
for radiative cooling, chemodynamics, and supernova feed- 
back. 

Our main conclusions are as follows: 

• The simulated stellar mass density of the spheroid 
(bulge+halo) falls off as r -3 ' 5 at large radii (r > 30 kpc), in 
agreement with the findings of previous theoretical studies 
based on hybrid analytic+N-body methods and cosmolog- 
ical hydrodynamical re-simulations, as well as deep obser- 
vations of M31. Fitting all ss 400 simulated galaxies sepa- 
rately, we find the 1-sigma (2-sigma) scatter in the logarith- 
mic slope is w 0.4 (~ 0.6). 

• In the inner region, the simulated mass density profile 
can be fit very well with a de Vaucouleurs profile with a 
median r e ff ~ 1.8 kpc (with considerable system-to-system 
scatter), and relatively well with a power-law with logarith- 
mic slope ~ — 3.0± —0.4. This is shallower than in the outer 
regions, but note that the mass density profile has a higher 
normalisation at small radii, so that a power-law with slope 
—3.5 extrapolated to the inner regions underpredicts the 
mass density there. The increased normalisation and still 
fairly steep profile at small radii differs from that found 
previously with hybrid analytic+N-body methods, which is 
most likely due to their neglect of in situ star formation. 
The shape of our simulated mass density profile matches 
that inferred for M31 (see Fig. [4} and other nearby galaxies 
remarkably well. 

• Significant negative metallicity gradients are a ubiqui- 
tous feature of the simulated galaxy populations. Typically 
[Fe/H] drops by 0.6 to 0.9 dex from the inner (r < 30 kpc) 
to the outer (r ~ 30 — 100 kpc) regions (see middle panel of 



Fig. [9} . T his drop is compatib le with that inferred for the 
MW (e.g-. lde Jong et al.ll2010l ) and M31 (see Fig.EJ. 

• On average, 2/3 of the mass of the stellar spheroid 
(bulge+halo) is composed of stars that formed in situ, but 
with considerable system-to-system scatter. Limiting the 
comparison to r > 20 kpc, this fraction drops to 1/5 on 
average. 

• The large contribution from in situ star formation in 
the inner regions drives the change in slope of the mass 
density profile at r ~ 30 kpc and also gives rise to large- 
scale metallicity gradients. 

• On average, the in situ stars were formed at radii com- 
parable to their present distance from the Galactic center, 
with considerable migration in both directions (towards and 
away from the center; see Fig. Hip . Typically, the stars were 
formed at z ~ 1 — 1.5, in a proto-disc. See McCarthy et al. 
(in prep.) for further discussion. 

• At fixed total mass, M200, galaxies with high D/T ratios 
actually have the most prominent stellar haloes (in terms 
of absolute mass density). This owes to the presence of a 
strong trend between D/T and total (disc+bulge+halo) stel- 
lar mass. 

In the Introduction we noted that the large differences 
in the properties of the stellar bulge+halo of the Milky Way 
and M31 (particularly their metallicities) has raised ques- 
tions about the compatibility with the ACDM scenario. As 
we have shown, however, there is considerable system-to- 
system scatter in the mean metallicities (as well as in the 
prominence of metallicity gradients) of the simulated Milky 
Way-mass systems (see Fig. . The scatter is roughly com- 
patible with the observed differences between the Milky Way 
and M31. 

M31 is typically consi dered to be more dynamica lly ac- 
tive than the Milky Way jHammer et alj|2007t l201Ch . The 
correlation between the epoch of formation and the in situ 
mass fraction found in the simulations would therefore sug- 
gest that the Milky Way has a higher in situ mas s frac- 
tion than M31 (as also argued bv lZolotov et al.ll2009f ). How- 
ever, this is difficult to reconcile with the fact that we find 
that the in situ mass fraction is also correlated with system 
metallicity, in the sense that more metal rich systems have 
higher in situ fractions (and that M31 has a more metal rich 
spheroidal component than the Milky Way) . However, as we 
have already noted, the simulations predict a large degree of 
scatter in these correlations, which may reconcile this issue. 

It would clearly be very interesting to measure the in 
situ mass fraction of the spheroidal components of the Milky 
Way and Andromeda and compare it directly to the predic- 
tions of the simulations (e.g., in Fig.[5J). Given their different 
formation histories, we expect there to be differences in the 
detailed spatial distributions and kinematics of the in situ 
and accreted stellar components, and this may hold the key 
to unambiguously separating these components from one an- 
other. We examine the kinematic properties, as well as the 
3D spatial structure, of the in situ and accreted components 
in McCarthy et al. (in prep.). 

It should be borne in mind that there are clear limi- 
tations in comparing a large sample of simulated galaxies 
(which could have unforseen "systematic errors", e.g., lack 
important physics) to one with very poor number statistics. 
Probing the diffuse stellar haloes of large samples of more 



Stellar haloes of disc galaxies 19 




Figure 12. Left: Median spherically-averaged stellar mass density profiles for the total bulge+halo component (solid) curves and just 
the in situ contribution to the bulge+halo component (dotted) curves. The black curves represent all 412 simulated galaxies, while the 
red curves represent only galaxies with 'hi' D/T (> 0.8) ratios (i.e., small fractional contributions from the bulge+halo) and the blue 
curves represent only galaxies with 'low' D/T (< 0.5) ratios (i.e., have a large fractional contribution from the bulge+halo). D/T has 
been computed using the method of CIO, described in Section 3.1. Right: Scatter plot of D/T versus total stellar mass. Systems with 
high D/T ratios actually have more prominent in situ stellar haloes in an absolute sense, since systems with high D/T ratios also tend 
to have higher total stellar masses. 



distant galaxies will go a very long way towards verifying 
(or refuting) the picture we have proposed. At present, it 
is only possible to probe the haloe s of distant galaxi es by 
stacking large samples of them (e.g. lZibetti et al.ll2004l ), but 
then one loses information about the system-to-system scat- 
ter, which is predicted to be large. In the coming years, ob- 
servations with the Large Synoptic Survey Telescope (for 
example), however, will allow us to probe the outer regions 
of distant galaxies in exquisite detail and will provide an 
excellent testbed for simulations of galaxy formation. 
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Figure 14. The mean satellite luminosity function of simu- 
lated galaxies in the intermediate-resolution (black) and high- 
resolution (red) GIMIC simulation. The hatched regions represent 
the Poisson error distributions for these luminosity function. For 
My < —12.5 there is good agreement. 

APPENDIX: NUMERICAL CONVERGENCE 

Here we investigate the sensitivity of our results to numerical 
resolution. As noted in Section 2, the —2a sphere has been 
simulated at eight times better mass resolution than the 
intermediate-resolution runs analysed in the present study. 

In Fig. [13] we compare the mass density and metallicity 
profiles (median + scatter) of 50 Milky Way-mass disc galax- 
ies in the high-resolution simulation with those of the ~ 400 
galaxies in the intermediate-resolution runs. Overall, the 
agreement between the high-resolution and intermediate- 
resolution GIMIC simulations is good for both the mass den- 
sity and metallicity profiles. 

We have also quantified the mass and fractional con- 
tribution of the in situ population to the spheroid in 
the high-resolution and intermediate-resolution simulations. 
The mean (median) mass of the in situ component is 
2.68 (2.48) xlO lo M and 2.01 (2.06) xl0 lo Af© in the 
intermediate- and high-res runs, respectively. Star forma- 
tion is therefore slightly less efficient in the high-res run 
(as can also been seen by the small difference in the nor- 
malisation of the density profiles in the top left panel of 
Fig. I13[) . The fractional contribution of the in situ popula- 
tion to the bulge+halo has a mean (median) of 0.71 (0.73) 
for the intermediate-resolution runs and has a mean (me- 
dian) of 0.69 (0.74) for the high-res run. 

In Fig. [14] we compare the mean satellite LF for 50 
Milky Way-mass disc galaxies in the high-resolution simu- 
lation with those of the « 400 galaxies in the intermediate- 
resolution runs. For My < —12.5 there is good convergence. 
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Figure 13. Left: Median spherically-averaged stellar mass density profiles for all Milky Way-mass systems in the high- and intermediate- 
resolution simulations. The scatter around the mean is indicated in each case by the different hashed regions (thinly spaced oblique lines 
for the intermediate-resolution run and widely spaced vertical lines for the high-resolution run). Right: Median spherically-averaged stellar 
metallicity profiles for the two runs. Overall, the agreement between the high-resolution and intermediate-resolution GIMIC simulations 
is good for both the mass density and metallicity profiles. 



